#include "base.h"
#include "acofs.h"

//carrega configuracoes do experimento
Experiment *newExperiment(Base *b, int antCount, int nElite, int b_local, unsigned long long int max, double minPher)
{
	Experiment *e;
	int i, j;

	e = (Experiment*) malloc(sizeof(Experiment));
	
	if(e)
	{
		e->elite = (char*) malloc(b->nFeatures * sizeof(char));
		e->nElite = nElite;
		
		e->bestSol = (char*) malloc(b->nFeatures * sizeof(char));
		e->bestSolF = -99999999;
		
		e->anthill = (Ant*) malloc(antCount * sizeof(Ant));
		e->nAnts = antCount;
		
		// Flag para realizar busca local
		e->b_local = b_local;
		
		e->trail = (double **) malloc(b->nFeatures * sizeof(double*));
		
		e->max = max;
		e->base = b;

		e->fInfoIndx = (int*) malloc(b->nFeatures * sizeof(int));
		indexFInfo(e->fInfoIndx, e->base->fInfo, b->nFeatures);

		if(!(e->elite && e->bestSol && e->anthill && e->trail &&  e->fInfoIndx)) {
			
			return NULL;
		}
		
		for(i = 0; i < b->nFeatures; i++)
		{
			e->elite[i] = 0;
			e->bestSol[i] = 0;
			
			e->trail[i] = (double *) malloc(b->nFeatures * sizeof(double));
            
			if(!e->trail[i])
				return NULL;
			
			for(j = 0; j < b->nFeatures; j++)
				e->trail[i][j] = (double) (abs(random() % 100)) / 1000 + minPher;
		}
		for(i = 0; i < e->nAnts; i++)
		{
            e->anthill[i].sol = (char *) malloc(b->nFeatures * sizeof(char));
			if(!e->anthill[i].sol) 
			{
				return NULL;
			}
		}

	}
	return e;
}

//limpa solucao da formiga
void resetAnt(Ant *a, int n)
{
	int i;
	
	if(a)
	{
		for(i = 0; i < n; i++) {
			
			a->sol[i] = 0;
		}
		
		a->solF = 0;
	}
}

//aco
void optimize(Experiment *e)
{

	long long unsigned int t = 0;
	int a;
	int eliteIndx[e->nAnts];
	long long unsigned int tBest = 0;
	long long unsigned int imprv = 0;

	int tmp, tmp1;
	while(t < e->max)
	{
		t = t + e->nAnts; //contabiliando antecipadamente
		for(a = 0; a < e->nAnts; a++)
		{
			//limpar formiga
			resetAnt(&(e->anthill[a]), e->base->nFeatures);
			
			//herdar da solucao elite
			inherit(e->anthill[a].sol, e->elite, e->base->nFeatures, 0.33);

			//construir caminho usando trail + heuristica construtiva
			buildway(&(e->anthill[a]), e->trail, e->base, e->fInfoIndx, 1.0, 1.0, 1.0);
			
			fprintf(stderr, "\nconstruido:\t");
			for(tmp = 0; tmp < e->base->nFeatures; tmp++)
				fprintf(stderr, "%d", e->anthill[a].sol[tmp] == 1 ? 1 : 0);
			
			fprintf(stderr, "\n");

			//avalia solução
			e->anthill[a].solF = h(e->base, 2, e->anthill[a].sol); // = funcao heuristica

			if(e->b_local) {
				
				t = localSearch(e->base, &(e->anthill[a]), t, e->max); //controle mais preciso de t na busca local
				fprintf(stderr, "\n buscado:\t");
				for(tmp = 0; tmp < e->base->nFeatures; tmp++)
					fprintf(stderr, "%d", e->anthill[a].sol[tmp] == 1 ? 1 : 0);
				fprintf(stderr, "\n");

			}
			
			fprintf(stderr, "\nformiga: %d|h: %f| inter %f | intra %f \n", a, e->anthill[a].solF, interClassDist(e->base, 2, e->anthill[a].sol), intraClassDist(e->base, 2, e->anthill[a].sol));
			
		}


    	//ordena solucoes
		rankAnts(e, eliteIndx);

		//eliminar solucao repetida?

		//seleciona elite
		elite(e, eliteIndx);

		//PRIMEIRA RODADA: usar pior formiga como referencia
		if(e->bestSolF <= -99999999)
		{
			e->bestSolF =  e->anthill[eliteIndx[e->nAnts - 1]].solF;
			saveBestSol(e, &(e->anthill[eliteIndx[e->nAnts - 1]]));
		}
		//atualiza feromonios
		updateTrail(e, eliteIndx);

		//guardar melhor solucao
        //fprintf(stderr, "melhor %f, melhor da rodada %f = formiga %d\n", e->bestSolF, e->anthill[eliteIndx[0]].solF, eliteIndx[0]);
        if(e->bestSolF < e->anthill[eliteIndx[0]].solF)
		{
			imprv++;
			tBest = t; //ultima melhora
            saveBestSol(e, &(e->anthill[eliteIndx[0]]));
		}

		fprintf(stderr, "trail\n");
		for(tmp = 0; tmp < e->base->nFeatures; tmp++)
		{
			for(tmp1 = 0;  tmp1 < e->base->nFeatures; tmp1++)
				fprintf(stderr, "%f ", e->trail[tmp][tmp1]);
			fprintf(stderr, "\n");
		}


    } //fim do while

	fprintf(stderr, "\ntBest %llu, imprv %llu\n", tBest, imprv);

}

//faz or das melhores features
void elite(Experiment *e, int *eliteIndx)
{
	int i, j;

	for(i = 0; i < e->base->nFeatures; i++)
	{
		e->elite[i] = 0;
		
		for(j = 0; j < e->nElite; j++) 
		{
			
			e->elite[i] = e->anthill[eliteIndx[j]].sol[i] || e->elite[i];
		}
	}
}

//ordena nao decrescente por Ant->solF indices das formigas
void rankAnts(Experiment *e, int *eliteIndx)
{
	int i, j, s, t;
	for(i = 0; i < e->nAnts; i++)
		eliteIndx[i] = i;
	
	for(i = 0; i < e->nAnts; i++)
	{
		s = i;
		for(j = i + 1; j < e->nAnts; j++)
		{
			if(e->anthill[eliteIndx[j]].solF > e->anthill[eliteIndx[s]].solF)
				s = j;
		}
		t = eliteIndx[i];
		eliteIndx[i] = eliteIndx[s];
		eliteIndx[s] = t;
	}
}

//gravar melhor solucao
void saveBestSol(Experiment *e, Ant *goldenAnt)
{
	int i;
	for(i = 0; i < e->base->nFeatures; i++)
		e->bestSol[i] = goldenAnt->sol[i];
	
	e->bestSolF = goldenAnt->solF;
}


//herdar solucao elite
//probabilidade p de herdar uma feature
void inherit(char *newSol, char *eliteSol, int nFeatures, double p)
{
	double r;
	int i;

	for(i = 0; i < nFeatures; i++)
		
		if(eliteSol[i])
		{
			r = (double) abs(random() % 1000000) / 1000000;
			
			if(r < p) 
			{
				newSol[i] = 1;
			}
		}
}

//evaporacao e adicao de feromonio
void updateTrail(Experiment *e, int *eliteIndx)
{
	int i, j;
	double v;
	for(i = 0; i < e->nElite; i++)
	{
		fprintf(stderr, "%f %f\n", e->bestSolF, e->anthill[eliteIndx[i]].solF);
		//valor de ajuste relacionado a posicao no ranking * qualidada da solucao em relacao a melhor solucao / constante de intensidade
		v = (e->nElite - i)  * (e->anthill[eliteIndx[i]].solF / e->bestSolF) / (e->nElite * 12);
		
		fprintf(stderr, "v = %f\n", v);
		for(j = 0; j < e->base->nFeatures; j++)
			if(e->anthill[eliteIndx[i]].sol[j] == 1)
				trailUp(j, e->trail, e->base->nFeatures, e->anthill[eliteIndx[i]].sol, v);

	}
	trailEvp(e->trail, e->base->nFeatures, 0.95);
}

//atualiza feromonio de quem esta junto de f na solucao
void trailUp(int f, double **trail, int n, char *sol, double v)
{
	int i;
	for(i = 0; i < n; i++)
		if(sol[i] == 1)
			trail[i][f] = trail[i][f] + v;

}

//evapora feromonio
void trailEvp(double **trail, int n, double v)
{
	int i, j;
	for(i = 0; i < n; i++)
		for(j = 0; j < n; j++)
		{
			trail[i][j] = trail[i][j] * v;
			if(trail[i][j] > 1)
				 trail[i][j] = 1;
			else if(trail[i][j] < 0)
				 trail[i][j] = 0;
		}
}

void indexFInfo(int *indx, double *v, int n)
{
	int i, j, k, t;

	for(i = 0; i < n; i++)
		indx[i] = i;

	for(i = 0; i < n; i++)
	{
		k = i;
		for(j = i + 1; j < n; j++)
			if(v[indx[j]] > v[indx[k]])
				k = j;
		t = indx[i];
		indx[i] = indx[k];
		indx[k] = t;
	}
}

//heuristica auxiliar na construcao de solucao
//usa info de feature e trilha
double hBuild(int f, int n, double **trail, double *fInfo, char *sol, double alpha, double beta)
{
	int i;
	int fs = 0;
	double v = 0;
	for(i = 0; i < n; i++)
		if(sol[i] == 1)
		{
			fs++;
			v = v + trail[i][f];
		}

	if(fs < 1)
		v = 0;
	else
		v = v / fs;
	
	return (alpha * v) + (beta * fInfo[f]);

}

//construcao de solucao
void buildway(Ant *a, double **trail, Base *b, int *indxFInfo, double alpha, double beta, double p) 
{
	int i, j;
	int fs = 0;
	int f = abs(random() % b->nFeatures);;
	double vf;
	double v = 0;
	
	//calcula v para features herdadas
	for(i = 0; i < b->nFeatures; i++)
		if(a->sol[i] == 1)
		{
			fs++;
			v = v + hBuild(i, b->nFeatures, trail, b->fInfo, a->sol, alpha, beta);
		}
	if(fs < 1)
		v = 0;
	else
		v = v / fs;

	//adcionar f?
	for(i = 0; i < b->nFeatures; i++)
	{
		f = (MAGIC_NUMBER + f ) % b->nFeatures; //supondo que nenhum problema tenha nFeatures multiplo de MAGIC_NUMBER
		//fprintf(stderr, "f %d == %d\n", f, a->sol[f]);
		if(a->sol[f] == 0) //ainda nao esta no conjunto
		{
			vf = hBuild(f, b->nFeatures, trail, b->fInfo, a->sol, alpha, beta);
			//fprintf(stderr, "vf %f - v %f \n", vf, v);
			if(vf >= v * p)
			{
			
				fs++;
				a->sol[f] = 1;
				trailUp(f, trail, b->nFeatures, a->sol, 0.007); //deposito de feromonio local
	
	
				v = 0; 
				//recalcular v para o conjunto atual de features
				for(j = 0; j < b->nFeatures; j++)
        			if(a->sol[j] == 1)
            			v = v + hBuild(j, b->nFeatures, trail, b->fInfo, a->sol, alpha, beta);
				
				v = v / fs;
			}
		}

	}
	trailEvp(trail, b->nFeatures, 0.99); //evaporacao local

}

double **calcClassCenters(Base *b, int numClasses, char *proj) {
	
	int i, j;
	int numElem[numClasses];
	long int class;
	
	// matriz contendo os centros de cada classe, as bases terao sempre 2 classes
	// positive e negative	
	double **center_vector = (double **) malloc(numClasses * sizeof(double *));
	if(!center_vector)
    {
        fprintf(stderr, "Malloc Fail::calcClassCenters\n");
        exit(-1);
    }
	
	for(i = 0; i < numClasses; i++) {
		
		numElem[i] = 0;
		center_vector[i] = (double *) malloc(b->nFeatures * sizeof(double));
		if(!center_vector[i])
    	{
        	fprintf(stderr, "Malloc Fail::calcClassCenters\n");
        	exit(-1);
    	}
		for(j = 0; j < b->nFeatures; j++)
			center_vector[i][j] = 0;

	}
	for(i = 0; i < b->nData; i++) 
	{
		class = (long int) b->database[i][b->nFeatures];
		numElem[class]++;
		for(j = 0; j < b->nFeatures; j++) 
		{
			if(proj[j] == 1)
				center_vector[class][j] += b->database[i][j];
		}
	}
	for (i = 0; i < numClasses; i++) 
	{
		for(j = 0; j < b->nFeatures; j++) 
		{
			if(proj[j] == 1)
				center_vector[i][j] = center_vector[i][j] / numElem[i];
			else
				center_vector[i][j] = 0;
		}
	}
	return center_vector;
}

double euclideanDist(double *vetA, double *vetB, int tamVet) {
	
	int i;
	double soma = 0.0;
	
	for (i = 0; i < tamVet; i ++) {
		
		soma += pow(vetA[i] - vetB[i], 2);
	}

	return sqrt(soma);
}

//Calcula distancia intra-classe
double interClassDist(Base *b, int numClasses, char *proj) 
{

	int i, j;
	double global_center_vector[b->nFeatures];
	double *result_dist = (double *) malloc (b->nData * sizeof(double));
	if(!result_dist)
	{
		fprintf(stderr, "Malloc Fail::interClassDist\n");
		exit(-1);
	}
	
	for(i = 0; i < b->nFeatures; global_center_vector[i++] = 0);

	//Calcula centro de cada classe
	double **center_vector = calcClassCenters(b, numClasses, proj);
	double d = 0;

	
	//Calcula o centro do conjunto de dados inteiro
	for(i = 0; i < b->nData; i++) 
	{					
		for(j = 0; j < b->nFeatures; j++) 
		{
			if(proj[j] == 1) 
				global_center_vector[j] += b->database[i][j]; //soma as colunas
		}
	}
	for(i = 0; i < b->nFeatures; i++) 
	{
			global_center_vector[i] = global_center_vector[i] / b->nData; //media das colunas (centro da base)
	}
	for (i = 0; i < numClasses; i++) {
	
		//result_dist[i] = 
		d += euclideanDist(global_center_vector, center_vector[i], b->nFeatures);
	}
	
	return d / numClasses;
}

double intraClassDist(Base *b, int numClasses, char *proj)
{
	int i;
	long int class;
	
	double *result_dist = (double *) malloc (b->nData * sizeof(double));	
	double **center_vector = calcClassCenters(b, numClasses, proj);
	double d = 0;
	for (i = 0; i < b->nData; i++) {
		
		class = (long int) b->database[i][b->nFeatures];
		//result_dist[i] = 
		d += euclideanDist((double *) b->database[i], center_vector[class], b->nFeatures);
	}
	return d / b->nData;
}

//funcao heuristica principal
double h(Base *b, int numClasses, char *proj)
{
	int i;
	int fs = 0;
	for(i = 0; i < b->nFeatures; i++)
		if(proj[i] == 1)
			fs++;
	//diminui o valor da funcao caso esteja selecionando mais que 80% das caracteristicas
	double v = 1;
	if(fs > b->nFeatures * 0.8)
		v = 1 + (0.8 - ((1.0 * fs) / (1.0 * b->nFeatures)));

	return v * ((1 * interClassDist(b, numClasses, proj)) / (1 *  intraClassDist(b, numClasses, proj)));
}

int localSearch(Base *b, Ant *a, int nEval, int MAX) 
{
	
	int i;
	int carac_selecionada = abs(random() % b->nFeatures);
	int fSelected;
	double v;
				
	while (nEval < MAX) {
		fSelected = 0;
		for(i = 0; i < b->nFeatures; i++)
			if(a->sol[i] == 1)
				fSelected++;
			
		//contribuicao media no fitness da solucao	
		v = a->solF / fSelected;

				
		// remove uma caracteristica
		i = 0;
		while(a->sol[carac_selecionada] != 1) {
		
			carac_selecionada = (carac_selecionada + MAGIC_NUMBER) % b->nFeatures;
			i++;
			if(i >= b->nFeatures) break;
		}
	
		a->sol[carac_selecionada] = 0;	
		double novo_valor_h = h(b, b->nClasses, a->sol);
		nEval++;
	
		fprintf(stderr, "fitness %f  \t novo fitness  %f\n", a->solF, novo_valor_h);
	
		if(novo_valor_h > a->solF) {

			fprintf(stderr, "removeu f%d\n", carac_selecionada);	
			a->solF = novo_valor_h;
			
			continue;
		}
		else {
		
			a->sol[carac_selecionada] = 1;
		}
	
		// adiciona uma caracteristica
		carac_selecionada = abs(random() % b->nFeatures);
		i = 0;
		
		i = 0;
		while(a->sol[carac_selecionada] != 0) {
		
			carac_selecionada = (carac_selecionada + MAGIC_NUMBER) % b->nFeatures;
			i++;
			if(i >= b->nFeatures) break;
		}
	
		a->sol[carac_selecionada] = 1;
		novo_valor_h = h(b, b->nClasses, a->sol);
		nEval++;
	
		fprintf(stderr, "fitness %f  \t novo fitness  %f\n", a->solF, novo_valor_h);
		if(novo_valor_h > a->solF) {
			fprintf(stderr, "adicionou  f%d\n", carac_selecionada);                                 
			a->solF = novo_valor_h;
			
			continue;
		}
		else {
		
			a->sol[carac_selecionada] = 0;
		}
	
		// Troca uma caracteristica aleatoria
		carac_selecionada = abs(random() % b->nFeatures);
	
		i = 0;
		while(a->sol[carac_selecionada] != 1) {
		
			carac_selecionada = (carac_selecionada + MAGIC_NUMBER) % b->nFeatures;
			i++;
			if(i >= b->nFeatures) break;
		}
	
		a->sol[carac_selecionada] = 0;
		int carac_selecionada2 = (carac_selecionada + MAGIC_NUMBER) % b->nFeatures;
		
		i = 0;
		while(a->sol[carac_selecionada2] != 0) {
		
			carac_selecionada2 = (carac_selecionada2 + MAGIC_NUMBER) % b->nFeatures;
			i++;
			if(i >= b->nFeatures) break;
		}
		
		a->sol[carac_selecionada2] = 1;

		novo_valor_h = h(b, b->nClasses, a->sol);
		nEval++;
	
		fprintf(stderr, "fitness %f  \t novo fitness  %f\n", a->solF, novo_valor_h);
		if(novo_valor_h > a->solF) {
		    fprintf(stderr, "removeu f%d e adicionou f%d\n", carac_selecionada, carac_selecionada2);                                 
			a->solF = novo_valor_h;
			
			continue;
		}
		else {
		
			a->sol[carac_selecionada] = 1;
			a->sol[carac_selecionada2] = 0;	
			return nEval;
		}
	}
}
